module jmcoef_df_mod; 
implicit none;
private;
public test_jmcoef;
contains;

subroutine test_jmcoef(); 
  integer:: ll   ! obital angular momentum
  real*8:: eta  ! k/lambda
  integer:: n   ! index n
  integer:: iz  ! charge
  real*8:: rsn, rcn   ! calculated coef.

  integer:: Nmax
  real*8:: theta, x, y

  ll = 0
  eta = 3.0/2.0
  Nmax = 40
  iz = 0

  theta = asin(eta/(eta*eta + 0.250d0))
  print*, 'theta=', theta
  x = (eta*eta - 0.250d0)/(eta*eta + 0.250d0)
  y = eta/(eta*eta + 0.250d0)
  print*, 'x,y=', x,y 

  do n=0,Nmax
     call jmcoef_df1(ll,eta,n,iz,rsn,1)

!todo     call jmcoef_df2(ll,eta,n,iz,rsn,rcn)
     print*, n, rsn, (-cos((n+1)*theta))/(n+1d0)!, (sin((n+1)*theta))/(n+1d0)
  enddo
end subroutine


!$$ This program calculate J-matrix expansion coefficeints for sine and cosine like solutions sn and cn.
!$$   D.Fursa, 17-04-2013
subroutine  jmcoef_df2(ll,rk,rlam,iz,n,rsn,rcn)
 
  USE HYP_2F1_MODULE

  implicit none
 
  integer, intent(in):: ll   ! orbital angular momentum
  real*8, intent(in):: rk  ! momentum
  real*8, intent(in):: rlam  ! Laguerre exp. fall-off
  integer, intent(in):: iz  ! Charge  (it is Z1 * Z2)
  integer, intent(in):: n   ! index n
  real*8, intent(out):: rsn, rcn   ! calculated coef., sin-like and cos-like

  real*8:: x,y, eta
  real*8:: tmp
  real*8 theta, t, Pi, coulphase
  complex(PR):: ztmp, ctmp1, ctmp2, ctmp, Ah, Bh, Ch, Zh, im, cphase_tmp
  complex(PR):: GAMMA_INV, LOG_GAMMA, HYP_2F1

  
  Pi = 3.141592654

  rsn = 0d0
  rcn = 0d0

  t = -dble(iz)/rk
  eta = rk/rlam

!!$  x = cos(theta); y = sin(theta)
  x = (eta*eta - 0.250d0)/(eta*eta + 0.250d0)
  y = eta/(eta*eta + 0.250d0)
  theta = asin(eta/(eta*eta + 0.250d0))
 

  im = (0,1d0)
  Ah = -ll - im*t
  Bh = n + 1d0
  Ch = n + ll + 2d0 - im*t
  Zh = exp(-2.0*im*theta)
  
  ctmp1 =  HYP_2F1(Ah,Bh,Ch,Zh)
  ztmp = n + ll + 2d0 - im*t
  ctmp2 = GAMMA_INV(ztmp)
  
  ztmp = ll + 1 - im*t
  cphase_tmp = exp(LOG_GAMMA(ztmp))
  coulphase = atan2(aimag(cphase_tmp),real(cphase_tmp))
  
  ztmp = n + 1
  tmp = exp(LOG_GAMMA(ztmp))
  
  ctmp = -ctmp1*ctmp2*tmp*exp(im*coulphase)*exp(Pi*t/2d0) * exp(-theta*t) * exp(-im*(n+1.0)*theta)/(2d0*y)**ll
  
  rcn = real(ctmp)   !   cos-like 
  rsn = aimag(ctmp)  !  sin-like
     
  return
end subroutine 


!// original by DmitryFursa 4-Apr-2013
!program test
!
!  implicit none 
!
!  integer:: ll   ! obital angular momentum
!  real*8:: eta  ! k/lambda
!  integer:: n   ! index n
!  integer:: iz  ! charge
!  real*8:: rsn, rcn   ! calculated coef.
!
!  integer:: Nmax
!  real*8:: theta, x, y
!
!  ll = 0
!  eta = 3.0/2.0
!  Nmax = 40
!  iz = 0
!
!  theta = asin(eta/(eta*eta + 0.250d0))
!  print*, 'theta=', theta
!  x = (eta*eta - 0.250d0)/(eta*eta + 0.250d0)
!  y = eta/(eta*eta + 0.250d0)
!  print*, 'x,y=', x,y 
!
!
!  do n=0,Nmax
!     call jmcoef(ll,eta,n,iz,rsn,1)
!     print*, n, rsn, (-cos((n+1)*theta))/(n+1d0)!, (sin((n+1)*theta))/(n+1d0)
!  enddo
!
!
!  stop 
!end program test


subroutine  jmcoef_df1(ll,eta,n,iz,rscn,iflag)
  implicit none 
  integer, intent(in):: ll   ! obital angular momentum
  real*8, intent(in):: eta  ! k/lambda
  integer, intent(in):: n   ! index n
  integer, intent(in):: iz  ! charge
  real*8, intent(out):: rscn   ! calculated coef.
  integer, intent(in):: iflag  ! =0 for sine like, =1 for cosine like 
  real*8:: x,y
  real*8:: u0, u1, un, unm1, unm2, ri, tmp, tt1, tt2, tt3, tt
  integer:: i, ittt
  !//integer:: igamma, iratiogamma, idoublefact
  !//real*8:: GausHyp

  
  rscn = 0d0
  
!!$  x = cos(theta); y = sin(theta)
  x = (eta*eta - 0.250d0)/(eta*eta + 0.250d0)
  y = eta/(eta*eta + 0.250d0)
!  print*, 'x,y=', x,y 

 
  if(iz .eq. 0) then
     
     if(iflag .eq. 0) then
        ittt = igamma(ll+1)
        u0 = 2**ll * ittt * y**(ll+1)
        u1 = 2*x*(ll+1) * u0
     elseif(iflag .eq. 1) then
        tt1 = idoublefact(ll)
        tt2 = igamma(2*ll+2)
        tt = (1d0-x)/2d0
        tt3 = GausHyp(ll,tt)    
        u0 = -(tt1/tt2)*(tt3/(sqrt(abs(tt)))**ll)
        u1 = 2*x*(ll+1) * u0 +idoublefact(ll+1)/y**ll
     else
        print*,' jmcoef(...): wrong value for iflag=', iflag
        stop
     endif
     print*, 'u0=', u0

     
     if(n .eq. 0) then

        un = u0

     elseif(n .eq. 1) then

        un = u1

     else

        unm2 = u0
        unm1 = u1
        do i=2,n
           ri = i
           un = (2*x*(i+ll)*unm1 - (i+2*ll)*unm2)/ri
!           print'(i5,3F15.5)',i, un, unm1, unm2
           unm2 = unm1
           unm1 = un
        enddo               

     endif

     tmp = iratiogamma(n,ll)
     rscn = un/tmp
     
  else
          
  endif
  return
end subroutine 

!$$----------------------------------

integer function igamma(n)
  implicit none
  integer, intent(in):: n
  integer:: j

  if(n .eq. 1) then
     igamma = 1
  elseif(n.eq. 2) then
     igamma = 1
  elseif(n.gt. 2) then
     
     igamma = 1
     do j=1,n-1
        igamma = igamma * j
     enddo
     
  else
     print*, 'igamma(...): Wrong argument: negative value or zero:',n
     stop
  endif


  return 
end function igamma
!$$----------------------------------

!$$ (2n-1)!! = \prod_{i=1}^{n} (2*n-1)
!$$
integer function idoublefact(n)
  implicit none
  integer, intent(in):: n
  integer:: j, itmp

     
  if(n .lt. 0) then
     print*, 'igamma(...): Wrong argument: negative value or zero:',n
     stop
  endif

  itmp = 1
  do j=1,n
     itmp = itmp * (2*j-1)
  enddo
     idoublefact = itmp

  return 
end function idoublefact

!$$----------------------------------

integer function iratiogamma(n,ll)
  implicit none
  integer, intent(in):: n, ll
  integer:: j, itmp

     
  if(n .lt. 0 .or. ll .lt. 0 ) then
     print*, 'iratiogamma(...): Wrong argument: negative values for n or ll:',n,ll
     stop
  endif

  itmp = 1
  do j=n+1,n+2*ll+1
     itmp = itmp * j
  enddo
     
  iratiogamma = itmp

  return 
end function iratiogamma

!$$----------------------------------
!$$ integer the Pochhammer symbol  = 1 if n=0; = a(a+1) ...(a+n-1)
integer function ipoch(i,n)
  implicit none
  integer, intent(in):: n, i
  integer:: j, itmp

     
  if(n .lt. 0 ) then
     print*, 'ipoch(...): Wrong argument: negative values for n:',n
     stop
  endif

  if( n .eq. 0) then
     ipoch = 1
  else
     ipoch = i
     do j=1,n-1
        ipoch = ipoch * (i + j)
     enddo     
  endif

  return 
end function ipoch

!$$ real the Pochhammer symbol  = 1 if n=0; = a(a+1) ...(a+n-1)
function ripoch(q,n)

  implicit none

  real*8:: ripoch
  real*8, intent(in):: q
  integer, intent(in):: n
  integer:: j, itmp

     
  if(n .lt. 0 ) then
     print*, 'ripoch(...): Wrong argument: negative values for n:',n
     stop
  endif

  if( n .eq. 0) then
     ripoch = 1
  else
     ripoch = q
     do j=1,n-1
        ripoch = ripoch * (q + j)
     enddo     
  endif

  return 
end function ripoch

!$$----------------------------------

!$$ This is Gauss hypergeometric Function:  2F1(-2*ll-1,1,1/2 - ll,tt)
!$$ In general 2F1(a,b,c,z) = \sum_{n=0}^{infty} \frac{(a)_n (b)_n}{(c)_n}  \frac{z^n}{n!}
!$$ where (a)_n is the Pochhammer symbol:  = 1 if n=0; = a(a+1) ...(a+n-1)  if n>0
!$$ Note that (1)_n = n!, or  ipoch(1,n) = igamma(n+1)

function GausHyp(ll,tt)

  implicit none
  
  real*8:: GausHyp
  integer, intent(in):: ll
  real*8, intent(in):: tt
  
  integer:: Nmax, n
  real*8:: tmp, ttn
  !real*8:: ripoch
  !integer:: ipoch

  tmp = 0
  Nmax = 2*ll + 1
  ttn = 1d0
  do n=0,Nmax
     tmp = tmp + (ipoch(-2*ll-1,n) / ripoch(0.5d0-ll,n))*ttn 
     ttn = ttn * tt
  enddo

  GausHyp = tmp

return
end function GausHyp

end module

